// fft.cu

// Copyright (C) 2010  traxex

// Author: 	traxex <traxex@yahoo.cn>
// Created: 	07 Jun 2010
// Version: 	1.0

#include "cuda/unittest.h"
#include <fftw3.h>

const int 		N 	= 16384;

void fftw_fft( cfloat *out, cfloat *in, int n){
    fftwf_plan p = fftwf_plan_dft_1d ( n, ( fftwf_complex* )in, ( fftwf_complex* )out, -1, FFTW_ESTIMATE);
    fftwf_execute ( p );    
    fftwf_destroy_plan ( p );
}

void fftw_ifft( cfloat *out, cfloat *in, int n){
    fftwf_plan p = fftwf_plan_dft_1d ( n, ( fftwf_complex* )in, ( fftwf_complex* )out, 1, FFTW_ESTIMATE);
    fftwf_execute ( p );    
    fftwf_destroy_plan ( p );
}


void performance( ){
    
    cfloat			h[N];
    vector<cfloat>		d(N);

    for( int i = 0; i < N; ++i)	h[i]=cfloat( i*i+i, i+2);
    d << h;

    //fftw ->complex[N]
    double t_fftw = _timeused({
            LOOP fftw_fft( h, h, N);
        });

    //cufft ->complex[N]
    double t_cufft = _timeused({
            LOOP fft( d, d);
        });
    printf( "cufft vs fftw:\t\t%.2f times\t(%d cfloat)\n", t_fftw/t_cufft, N);

    matrix<cfloat>	d2(LOOPN,N);
    double t_nfft = _timeused({
            fft( d2, d2);
        });
    printf( "cufft many vs cufft:\t%.2f times\t(%d*%d cfloat)\n", t_cufft/t_nfft, LOOPN, N);
}


bool unittest(){
    const int N = 20;
    cfloat in[N];
    
    cfloat out1[N];
    cfloat out2[N];

    vector<cfloat> din(N);
    vector<cfloat> dout(N);

    //fftw init
    fftw_fft( out1, out1, N);
    fftw_ifft( out2, out2, N);
    
    //init 
    for( int i = 0; i < N; ++i){ in[i] = cfloat(i*i+i, 3*i+10); }
    din << in;

    fftw_fft( out1, in, N);					//fftw fft
    fft( dout, din);						//cuda fft
    dout >> out2;

    //check
    if( !vcmp( out1, out2, N)) return false;

    fftw_ifft( out1, in, N);    				//fftw ifft
    ifft( dout, din);						//cuda ifft
    dout >> out2;

    if( !vcmp( out1, out2, N)) return false;
    
    return true;
}


